Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Origin

8.3.1.1 Positively associated genomes with domestic cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_CaboVerde.50	Actinomycetota	Actinomycetia	Actinomycetales	Actinomycetaceae	NA	0.953
bin_Aruba.25	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium pseudocatenulatum	0.995
bin_Aruba.2	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium adolescentis	0.994
bin_Aruba.6	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium longum	0.985
bin_Denmark.14	Actinomycetota	Actinomycetia	Actinomycetales	Bifidobacteriaceae	Bifidobacterium gallinarum	0.950
bin_Aruba.47	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	NA	0.964
bin_Aruba.57	Actinomycetota	Actinomycetia	Mycobacteriales	Mycobacteriaceae	Corynebacterium pyruviciproducens	0.931
bin_Aruba.51	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	Parolsenella sp900544655	0.969
bin_Denmark.44	Actinomycetota	Coriobacteriia	Coriobacteriales	Atopobiaceae	UBA7748 sp900314535	0.951
bin_Aruba.14	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Denmark.4	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	1.000
bin_Aruba.21	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	NA	0.999
bin_Brazil.61	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella tanakaei	0.998
bin_Aruba.39	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp902362275	0.992
bin_Spain.84	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900555555	0.992
bin_Brazil.151	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella intestinalis	0.991
bin_Denmark.127	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella sp900548365	0.991
bin_Denmark.17	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella bouchesdurhonensis	0.983
bin_Brazil.81	Bacillota	Bacilli	RFN20	CAG-826	NA	0.907
bin_Brazil.109	Bacillota_A	Clostridia	Oscillospirales	Butyricicoccaceae	AM07-15 sp003477405	0.936
bin_Aruba.28	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.995
bin_Malaysia.9	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Dysosmobacter welbionis	0.986
bin_Malaysia.116	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Flavonifractor plautii	0.981
bin_Malaysia.34	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Lawsonibacter sp000177015	0.954
bin_Aruba.54	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.959
bin_CaboVerde.35	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.930
bin_Aruba.36	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	Peptoniphilus_C sp902363535	0.913
bin_Aruba.52	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.906
bin_Malaysia.26	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	NA	0.976
bin_Aruba.29	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Gemmiger variabilis_C	0.941
bin_Malaysia.103	Bacillota_C	Negativicutes	Acidaminococcales	Acidaminococcaceae	Acidaminococcus intestini	0.942
bin_Malaysia.8	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Dialister sp900543165	0.981
bin_Denmark.60	Bacillota_C	Negativicutes	Veillonellales	Dialisteraceae	Allisonella histaminiformans	0.938
bin_Aruba.64	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.970
bin_Brazil.58	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera elsdenii	0.969
bin_CaboVerde.38	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera sp000417505	0.961
bin_Malaysia.81	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Megasphaera stantonii	0.943
bin_Malaysia.96	Bacillota_C	Negativicutes	Veillonellales	Megasphaeraceae	Caecibacter sp003467125	0.919
bin_Brazil.163	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella lascolaii	0.992
bin_CaboVerde.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella copri	0.981
bin_Aruba.10	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900544825	0.978
bin_Brazil.5	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900540415	0.929
bin_Spain.21	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprophilus	0.929
bin_Malaysia.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp000437675	0.924
bin_Malaysia.117	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp900541335	0.919
bin_Malaysia.77	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900542985	0.913
bin_Brazil.75	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.964
bin_Malaysia.61	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	Helicobacter_C labetoulli	0.964
bin_Brazil.128	Campylobacterota	Campylobacteria	Campylobacterales	Helicobacteraceae	NA	0.963
bin_Brazil.70	Desulfobacterota	Desulfovibrionia	Desulfovibrionales	Desulfovibrionaceae	Mailhella sp003150275	0.944
bin_Brazil.146	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp001917125	0.969
bin_Brazil.9	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp000436375	0.919
bin_Aruba.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.963
bin_Brazil.51	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	Sutterella wadsworthensis_A	0.963
bin_Spain.43	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.959
bin_Brazil.64	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.948
bin_Brazil.92	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	CAG-521 sp000437635	0.948
bin_Brazil.15	Pseudomonadota	Gammaproteobacteria	Burkholderiales	Burkholderiaceae	NA	0.930
bin_CaboVerde.33	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum sp900543125	0.995
bin_Brazil.82	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum succiniciproducens	0.990
bin_CaboVerde.55	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum_A thomasii	0.939

8.3.1.2 Positively associated genomes feral cats

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_Denmark.27	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_B hirae	0.000
bin_Brazil.12	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus faecalis	0.001
bin_Brazil.170	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	NA	0.014
bin_CaboVerde.24	Bacillota	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus_E cecorum	0.020
bin_CaboVerde.16	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_Denmark.56	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Levilactobacillus brevis	0.000
bin_Malaysia.72	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Limosilactobacillus reuteri	0.000
bin_CaboVerde.60	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Malaysia.127	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus agilis	0.001
bin_Denmark.61	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Latilactobacillus sakei	0.003
bin_Malaysia.35	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus animalis	0.003
bin_Malaysia.4	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus ruminis	0.035
bin_CaboVerde.14	Bacillota	Bacilli	Lactobacillales	Lactobacillaceae	Lactobacillus acidophilus	0.067
bin_Aruba.18	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus lutetiensis	0.001
bin_Brazil.22	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus pasteurianus	0.001
bin_Denmark.117	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Lactococcus garvieae	0.001
bin_Denmark.113	Bacillota	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus alactolyticus	0.002
bin_Brazil.69	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Clostridium_A leptum	0.097
bin_Denmark.93	Bacillota_A	Clostridia	Lachnospirales	Anaerotignaceae	Anaerotignum sp001304995	0.074
bin_CaboVerde.3	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Denmark.42	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium_P perfringens	0.000
bin_Brazil.136	Bacillota_A	Clostridia	Clostridiales	Clostridiaceae	Clostridium thermobutyricum	0.041
bin_Brazil.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900555395	0.026
bin_Brazil.89	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900550235	0.034
bin_Denmark.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerostipes hadrus	0.040
bin_Denmark.19	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900556835	0.052
bin_Denmark.118	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_A longicatena	0.067
bin_Malaysia.108	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.081
bin_Brazil.113	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerobutyricum sp900754855	0.095
bin_Denmark.50	Bacillota_A	Clostridia	Peptostreptococcales	Peptostreptococcaceae	Peptacetobacter sp900539645	0.053
bin_Brazil.107	Bacillota_A	Clostridia	Monoglobales_A	UBA1381	CAG-41 sp900066215	0.016
bin_Brazil.159	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900541465	0.040
bin_Brazil.140	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_A sp900543175	0.046
bin_Malaysia.6	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900545035	0.088

8.3.1.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.1.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)

***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.3.1.5 Predict responses (origin)

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")

# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.3.1.6 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.3.1.6.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0103	0.008668973	0.0005965939	0.016001173	0.915	0.085	Biosynthesis	Nucleic acid biosynthesis	UDP/UTP
D0504	0.004647085	0.0003771070	0.009768176	0.922	0.078	Degradation	Amino acid degradation	Methionine
D0906	0.003822519	0.0002160812	0.008201395	0.932	0.068	Degradation	Antibiotic degradation	Fosfomycin
D0205	0.012204272	0.0019631234	0.022408114	0.950	0.050	Degradation	Polysaccharide degradation	Pectin
D0208	0.009799367	0.0014936231	0.017659224	0.939	0.061	Degradation	Polysaccharide degradation	Mixed-Linkage glucans

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.004132651	-0.009242449	-2.704579e-04	0.044	0.956	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021408344	-0.038627047	-4.942490e-03	0.059	0.941	Biosynthesis	Amino acid biosynthesis	Glutamate
B0204	-0.016044178	-0.032154456	-1.338553e-03	0.082	0.918	Biosynthesis	Amino acid biosynthesis	Serine
B0302	-0.004993799	-0.010547888	-6.852718e-04	0.032	0.968	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.012631605	-0.022578453	-3.174742e-03	0.035	0.965	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011515181	-0.020902183	-2.246292e-03	0.066	0.934	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0309	-0.007441518	-0.015073476	-1.638992e-04	0.095	0.905	Biosynthesis	Amino acid derivative biosynthesis	Putrescine
B0804	-0.016343376	-0.029354205	-3.591691e-03	0.049	0.951	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016580393	-0.031741463	-2.636490e-03	0.060	0.940	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009240453	-0.017544737	-1.391204e-03	0.069	0.931	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.011102735	-0.022069357	-2.102256e-04	0.094	0.906	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002073624	-0.003638992	-6.474296e-04	0.033	0.967	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004495743	-0.007874479	-1.353705e-03	0.023	0.977	Degradation	Amino acid degradation	Ornithine
D0508	-0.003245058	-0.007138266	-1.816237e-04	0.079	0.921	Degradation	Amino acid degradation	Lysine
D0512	-0.018089246	-0.035972853	-1.186182e-03	0.086	0.914	Degradation	Amino acid degradation	Histidine
D0903	-0.004117914	-0.009256580	-2.622164e-04	0.044	0.956	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015447072	-0.027438899	-3.764076e-03	0.056	0.944	Degradation	Antibiotic degradation	Macrolide
D0601	-0.009517923	-0.017285585	-2.651368e-03	0.024	0.976	Degradation	Nitrogen compound degradation	Nitrate
D0603	-0.001981904	-0.003844551	-3.973787e-04	0.036	0.964	Degradation	Nitrogen compound degradation	Urate
D0610	-0.002955038	-0.004863991	-1.021881e-03	0.039	0.961	Degradation	Nitrogen compound degradation	Methylamine
D0611	-0.004117914	-0.009256580	-2.622164e-04	0.044	0.956	Degradation	Nitrogen compound degradation	Phenylethylamine
D0606	-0.005881451	-0.011632320	-9.896364e-04	0.063	0.937	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001574312	-0.002955348	-1.518534e-04	0.076	0.924	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001305868	-0.002159595	-9.074796e-05	0.008	0.992	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001305868	-0.002159595	-9.074796e-05	0.008	0.992	Degradation	Xenobiotic degradation	Xylene
D0807	-0.004098974	-0.008643808	-7.646272e-04	0.043	0.957	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004953919	-0.010591016	-6.720724e-04	0.049	0.951	Degradation	Xenobiotic degradation	Trans-cinnamate
D0816	-0.005871819	-0.011500014	-4.581529e-04	0.083	0.917	Degradation	Xenobiotic degradation	Phenylacetate

8.3.1.6.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()

trait	mean	p1	p9	positive_support	negative_support	Domain	Function	Element
B0219	-0.004132651	-0.009242449	-2.704579e-04	0.044	0.956	Biosynthesis	Amino acid biosynthesis	GABA
B0214	-0.021408344	-0.038627047	-4.942490e-03	0.059	0.941	Biosynthesis	Amino acid biosynthesis	Glutamate
B0204	-0.016044178	-0.032154456	-1.338553e-03	0.082	0.918	Biosynthesis	Amino acid biosynthesis	Serine
B0302	-0.004993799	-0.010547888	-6.852718e-04	0.032	0.968	Biosynthesis	Amino acid derivative biosynthesis	Betaine
B0310	-0.012631605	-0.022578453	-3.174742e-03	0.035	0.965	Biosynthesis	Amino acid derivative biosynthesis	Tryptamine
B0303	-0.011515181	-0.020902183	-2.246292e-03	0.066	0.934	Biosynthesis	Amino acid derivative biosynthesis	Ectoine
B0309	-0.007441518	-0.015073476	-1.638992e-04	0.095	0.905	Biosynthesis	Amino acid derivative biosynthesis	Putrescine
B0804	-0.016343376	-0.029354205	-3.591691e-03	0.049	0.951	Biosynthesis	Aromatic compound biosynthesis	Dipicolinate
B0603	-0.016580393	-0.031741463	-2.636490e-03	0.060	0.940	Biosynthesis	Organic anion biosynthesis	Citrate
B0601	-0.009240453	-0.017544737	-1.391204e-03	0.069	0.931	Biosynthesis	Organic anion biosynthesis	Succinate
B0401	-0.011102735	-0.022069357	-2.102256e-04	0.094	0.906	Biosynthesis	SCFA biosynthesis	Acetate
B0709	-0.002073624	-0.003638992	-6.474296e-04	0.033	0.967	Biosynthesis	Vitamin biosynthesis	Tocopherol/tocotorienol (E)
D0517	-0.004495743	-0.007874479	-1.353705e-03	0.023	0.977	Degradation	Amino acid degradation	Ornithine
D0508	-0.003245058	-0.007138266	-1.816237e-04	0.079	0.921	Degradation	Amino acid degradation	Lysine
D0512	-0.018089246	-0.035972853	-1.186182e-03	0.086	0.914	Degradation	Amino acid degradation	Histidine
D0903	-0.004117914	-0.009256580	-2.622164e-04	0.044	0.956	Degradation	Antibiotic degradation	Cephalosporin
D0908	-0.015447072	-0.027438899	-3.764076e-03	0.056	0.944	Degradation	Antibiotic degradation	Macrolide
D0601	-0.009517923	-0.017285585	-2.651368e-03	0.024	0.976	Degradation	Nitrogen compound degradation	Nitrate
D0603	-0.001981904	-0.003844551	-3.973787e-04	0.036	0.964	Degradation	Nitrogen compound degradation	Urate
D0610	-0.002955038	-0.004863991	-1.021881e-03	0.039	0.961	Degradation	Nitrogen compound degradation	Methylamine
D0611	-0.004117914	-0.009256580	-2.622164e-04	0.044	0.956	Degradation	Nitrogen compound degradation	Phenylethylamine
D0606	-0.005881451	-0.011632320	-9.896364e-04	0.063	0.937	Degradation	Nitrogen compound degradation	Allantoin
D0612	-0.001574312	-0.002955348	-1.518534e-04	0.076	0.924	Degradation	Nitrogen compound degradation	Hypotaurine
D0801	-0.001305868	-0.002159595	-9.074796e-05	0.008	0.992	Degradation	Xenobiotic degradation	Toluene
D0802	-0.001305868	-0.002159595	-9.074796e-05	0.008	0.992	Degradation	Xenobiotic degradation	Xylene
D0807	-0.004098974	-0.008643808	-7.646272e-04	0.043	0.957	Degradation	Xenobiotic degradation	Catechol
D0817	-0.004953919	-0.010591016	-6.720724e-04	0.049	0.951	Degradation	Xenobiotic degradation	Trans-cinnamate
D0816	-0.005871819	-0.011500014	-4.581529e-04	0.083	0.917	Degradation	Xenobiotic degradation	Phenylacetate

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.1.7 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D02	8.252038e-03	-0.0029397646	0.0196760509	0.822	0.178
B08	7.109220e-03	-0.0040068882	0.0171378858	0.770	0.230
B01	1.216320e-03	-0.0059241188	0.0081540077	0.610	0.390
S01	1.014944e-03	-0.0109145072	0.0133945495	0.551	0.449
B09	8.119340e-07	-0.0005326007	0.0004986745	0.367	0.633

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D08	-1.093511e-03	-0.002214746	-2.281340e-04	0.032	0.968
B03	-1.057672e-02	-0.017631148	-3.229192e-03	0.045	0.955
D06	-3.229620e-03	-0.006962169	-5.139396e-05	0.098	0.902
B04	-7.665435e-03	-0.017979073	1.779859e-03	0.149	0.851
B06	-6.778849e-03	-0.016584138	2.083246e-03	0.171	0.829
D07	-1.166506e-02	-0.029621342	4.098672e-03	0.194	0.806
D03	-4.540739e-03	-0.012618543	2.717105e-03	0.207	0.793
D05	-1.909669e-03	-0.007386594	2.833109e-03	0.221	0.779
B02	-3.789939e-03	-0.012493886	4.029648e-03	0.247	0.753
S03	-9.214706e-03	-0.031916056	1.659255e-02	0.249	0.751
D09	-1.644702e-03	-0.007831696	4.607589e-03	0.318	0.682
S02	-4.283466e-03	-0.013942829	2.781315e-03	0.319	0.681
B07	-3.643201e-03	-0.015378159	8.373293e-03	0.322	0.678
D01	-3.039883e-04	-0.004674591	4.160412e-03	0.417	0.583
B10	-1.661264e-05	-0.000312485	2.497855e-04	0.449	0.551

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2 Sex

8.3.2.1 Negatively associated genomes with male cats

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()

variable	mean	sd
Random: location	37.900015	25.317903
logseqdepth	56.110626	25.796874
sex	4.937460	5.612719
origin	1.051899	1.282563

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    column_to_rownames(var="genome")

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="sexMale", trend=="Negative") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()

genome	phylum	class	order	family	species	value
bin_Denmark.4	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	0.090
bin_Aruba.14	Actinomycetota	Coriobacteriia	Coriobacteriales	Coriobacteriaceae	Collinsella stercoris	0.089
bin_Brazil.91	Bacillota	Bacilli	RF39	UBA660	CAG-988 sp003149915	0.082
bin_Brazil.76	Bacillota	Bacilli	RF39	UBA660	CAG-877 sp900554305	0.017
bin_Malaysia.17	Bacillota	Bacilli	RF39	UBA660	NA	0.016
bin_Brazil.45	Bacillota	Bacilli	RF39	UBA660	CAG-877 sp900554315	0.010
bin_Malaysia.16	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	NA	0.085
bin_Malaysia.78	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Ruminococcus_E bromii_B	0.081
bin_Brazil.69	Bacillota_A	Clostridia	Oscillospirales	Acutalibacteraceae	Clostridium_A leptum	0.020
bin_Brazil.83	Bacillota_A	Clostridia	Lachnospirales	Anaerotignaceae	NA	0.069
bin_Denmark.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerostipes hadrus	0.091
bin_Spain.11	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.084
bin_Malaysia.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.077
bin_Malaysia.45	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A wexlerae	0.077
bin_Brazil.1	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.071
bin_Malaysia.7	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	CAG-317 sp900543415	0.070
bin_Denmark.118	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_A longicatena	0.066
bin_Brazil.3	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900555395	0.065
bin_Malaysia.52	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Fusicatenibacter saccharivorans	0.065
bin_Brazil.89	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900550235	0.062
bin_Denmark.34	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Fusicatenibacter saccharivorans	0.062
bin_Malaysia.73	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.062
bin_Brazil.166	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.057
bin_Brazil.165	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A caecimuris	0.055
bin_Denmark.19	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Faecalimonas sp900556835	0.053
bin_Spain.53	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.053
bin_Brazil.54	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Dorea_B phocaeensis	0.051
bin_Brazil.157	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.046
bin_Brazil.32	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Mediterraneibacter torques	0.042
bin_Brazil.8	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia_A sp900547615	0.042
bin_Spain.6	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Schaedlerella glycyrrhizinilytica	0.039
bin_Denmark.109	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Eisenbergiella sp900550285	0.034
bin_Brazil.113	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Anaerobutyricum sp900754855	0.031
bin_Brazil.17	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.030
bin_Spain.37	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	UMGS1472 sp900552095	0.029
bin_Denmark.52	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_B gnavus	0.023
bin_Brazil.56	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_A sp000432335	0.022
bin_Spain.67	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Ruminococcus_B sp900544395	0.021
bin_Brazil.28	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.019
bin_Brazil.93	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.019
bin_Malaysia.110	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Catenibacillus sp900557175	0.018
bin_Aruba.43	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Roseburia sp900548205	0.017
bin_Malaysia.108	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.017
bin_Brazil.116	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.016
bin_Malaysia.86	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.015
bin_Brazil.99	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospira sp900552795	0.013
bin_CaboVerde.61	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Roseburia sp900548205	0.013
bin_Brazil.97	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	UBA9502 sp900538475	0.010
bin_Malaysia.46	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp003287895	0.010
bin_Brazil.50	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospira sp003451515	0.009
bin_Denmark.66	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp900120295	0.009
bin_Spain.60	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.009
bin_Malaysia.125	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Blautia sp900539145	0.008
bin_CaboVerde.34	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.005
bin_Aruba.13	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	CAG-81 sp000435795	0.002
bin_Brazil.63	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.002
bin_Denmark.62	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Clostridium_Q sp000435655	0.002
bin_Brazil.105	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	NA	0.001
bin_Denmark.20	Bacillota_A	Clostridia	Lachnospirales	Lachnospiraceae	Enterocloster sp001517625	0.001
bin_Malaysia.98	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.066
bin_Malaysia.9	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Dysosmobacter welbionis	0.065
bin_Aruba.28	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.044
bin_Brazil.177	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.044
bin_Malaysia.116	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Flavonifractor plautii	0.043
bin_Malaysia.34	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	Lawsonibacter sp000177015	0.040
bin_Brazil.137	Bacillota_A	Clostridia	Oscillospirales	Oscillospiraceae	NA	0.027
bin_Aruba.52	Bacillota_A	Clostridia	Tissierellales	Peptoniphilaceae	NA	0.094
bin_Aruba.31	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Negativibacillus sp000435195	0.080
bin_Malaysia.102	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	NA	0.080
bin_Malaysia.30	Bacillota_A	Clostridia	Oscillospirales	Ruminococcaceae	Phocea massiliensis	0.055
bin_Denmark.72	Bacillota_C	Negativicutes	Selenomonadales	Selenomonadaceae	Megamonas funiformis	0.082
bin_Denmark.85	Bacillota_C	Negativicutes	Selenomonadales	Selenomonadaceae	NA	0.077
bin_Brazil.5	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900540415	0.080
bin_Malaysia.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp000437675	0.051
bin_Brazil.48	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides sp900766005	0.050
bin_Malaysia.117	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotellamassilia sp900541335	0.047
bin_CaboVerde.18	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella copri	0.038
bin_Aruba.10	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella sp900544825	0.037
bin_Malaysia.64	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Paraprevotella clara	0.033
bin_Brazil.163	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Prevotella lascolaii	0.029
bin_Spain.4	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp000436795	0.027
bin_Spain.48	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides fragilis	0.019
bin_Brazil.96	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides uniformis	0.014
bin_Malaysia.105	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola massiliensis	0.014
bin_Spain.21	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprophilus	0.012
bin_Malaysia.121	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides stercoris	0.011
bin_Brazil.43	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola vulgatus	0.010
bin_Denmark.30	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900546645	0.010
bin_Brazil.103	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola plebeius_A	0.009
bin_Malaysia.77	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola sp900542985	0.006
bin_Brazil.6	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Phocaeicola coprocola	0.005
bin_Brazil.110	Bacteroidota	Bacteroidia	Bacteroidales	Barnesiellaceae	Barnesiella intestinihominis	0.024
bin_Malaysia.13	Bacteroidota	Bacteroidia	Bacteroidales	Muribaculaceae	NA	0.075
bin_Spain.23	Bacteroidota	Bacteroidia	Bacteroidales	Muribaculaceae	CAG-279 sp900541935	0.051
bin_CaboVerde.37	Bacteroidota	Bacteroidia	Bacteroidales	Porphyromonadaceae	NA	0.058
bin_Brazil.11	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes putredinis	0.045
bin_Malaysia.131	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes putredinis	0.043
bin_Brazil.38	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes communis	0.039
bin_Brazil.86	Bacteroidota	Bacteroidia	Bacteroidales	Rikenellaceae	Alistipes dispar	0.030
bin_Brazil.138	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	Parabacteroides sp000436495	0.023
bin_Brazil.124	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	NA	0.019
bin_Brazil.160	Bacteroidota	Bacteroidia	Bacteroidales	Tannerellaceae	Parabacteroides distasonis	0.013
bin_CaboVerde.10	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	NA	0.070
bin_Brazil.68	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	Campylobacter_D helveticus	0.049
bin_Spain.26	Campylobacterota	Campylobacteria	Campylobacterales	Campylobacteraceae	Campylobacter_D upsaliensis	0.047
bin_Brazil.145	Cyanobacteria	Vampirovibrionia	Gastranaerophilales	Gastranaerophilaceae	NA	0.044
bin_Brazil.140	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_A sp900543175	0.029
bin_Malaysia.6	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900545035	0.015
bin_Brazil.159	Fusobacteriota	Fusobacteriia	Fusobacteriales	Fusobacteriaceae	Fusobacterium_B sp900541465	0.012
bin_Brazil.146	Pseudomonadota	Alphaproteobacteria	RF32	CAG-239	CAG-495 sp001917125	0.074
bin_Malaysia.50	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	NA	0.084
bin_Brazil.82	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum succiniciproducens	0.046
bin_CaboVerde.33	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum sp900543125	0.045
bin_CaboVerde.55	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Anaerobiospirillum_A thomasii	0.036
bin_Brazil.111	Pseudomonadota	Gammaproteobacteria	Enterobacterales	Succinivibrionaceae	Succinivibrio sp000431835	0.024

8.3.2.2 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="sexMale") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()

inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.3.2.3 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("Male","Female","Unknown")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "sex", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(sex=rep(gradient,1000)) %>%
            pivot_longer(!sex,names_to = "genome", values_to = "value")

# weights:  4 (3 variable)
initial  value 63.769541 
final  value 61.728471 
converged

8.3.2.4 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.3.2.4.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)

# A tibble: 30 × 9
   trait    mean       p1      p9 positive_support negative_support Domain       Function                           Element       
   <chr>   <dbl>    <dbl>   <dbl>            <dbl>            <dbl> <chr>        <chr>                              <chr>         
 1 B0220 0.00893 0.00137  0.0165             0.924            0.076 Biosynthesis Amino acid biosynthesis            Beta-alanine  
 2 B0212 0.0293  0.00379  0.0522             0.922            0.078 Biosynthesis Amino acid biosynthesis            Arginine      
 3 B0218 0.0110  0.000319 0.0224             0.902            0.098 Biosynthesis Amino acid biosynthesis            Tyrosine      
 4 B0310 0.0183  0.00536  0.0312             0.964            0.036 Biosynthesis Amino acid derivative biosynthesis Tryptamine    
 5 B0303 0.0145  0.00204  0.0255             0.924            0.076 Biosynthesis Amino acid derivative biosynthesis Ectoine       
 6 B0307 0.0418  0.00475  0.0711             0.914            0.086 Biosynthesis Amino acid derivative biosynthesis Spermidine    
 7 B0901 0.00118 0.000115 0.00255            0.917            0.083 Biosynthesis Metallophore biosynthesis          Staphyloferrin
 8 B0105 0.0262  0.0106   0.0428             0.948            0.052 Biosynthesis Nucleic acid biosynthesis          ADP/ATP       
 9 B0104 0.0351  0.0102   0.0581             0.942            0.058 Biosynthesis Nucleic acid biosynthesis          CDP/CTP       
10 B0106 0.0150  0.00305  0.0267             0.921            0.079 Biosynthesis Nucleic acid biosynthesis          GDP/GTP       
# ℹ 20 more rows

8.3.2.4.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)

# A tibble: 22 × 9
   trait     mean       p1        p9 positive_support negative_support Domain       Function                   Element           
   <chr>    <dbl>    <dbl>     <dbl>            <dbl>            <dbl> <chr>        <chr>                      <chr>             
 1 B0216 -0.0486  -0.0729  -0.0249              0.026            0.974 Biosynthesis Amino acid biosynthesis    Tryptophan        
 2 B0215 -0.0591  -0.0903  -0.0279              0.03             0.97  Biosynthesis Amino acid biosynthesis    Histidine         
 3 B0208 -0.0292  -0.0541  -0.00505             0.069            0.931 Biosynthesis Amino acid biosynthesis    Valine            
 4 B0209 -0.0283  -0.0555  -0.00322             0.073            0.927 Biosynthesis Amino acid biosynthesis    Isoleucine        
 5 B1012 -0.00833 -0.0130  -0.00368             0.005            0.995 Biosynthesis Antibiotic biosynthesis    Fosfomycin        
 6 B1004 -0.00460 -0.00875 -0.000640            0.081            0.919 Biosynthesis Antibiotic biosynthesis    Bacilysin         
 7 B0604 -0.0259  -0.0515  -0.000348            0.098            0.902 Biosynthesis Organic anion biosynthesis L-lactate         
 8 B0704 -0.0583  -0.0890  -0.0273              0.013            0.987 Biosynthesis Vitamin biosynthesis       Pantothenate (B5) 
 9 B0711 -0.0330  -0.0522  -0.0149              0.041            0.959 Biosynthesis Vitamin biosynthesis       Menaquinone (K2)  
10 B0710 -0.0179  -0.0297  -0.00596             0.057            0.943 Biosynthesis Vitamin biosynthesis       Phylloquinone (K1)
# ℹ 12 more rows

positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
#      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.3.2.5 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(sex, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "sex") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="sex")
   })

#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "sex") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
D07	0.0375944479	1.972014e-02	0.056148600	0.992	0.008
B01	0.0120235113	2.924545e-03	0.021231683	0.939	0.061
B03	0.0134706279	3.129619e-03	0.021856061	0.932	0.068
D06	0.0046369186	1.904574e-04	0.009455340	0.909	0.091
B09	0.0008808116	2.773819e-05	0.001863172	0.905	0.095
D03	0.0086358818	-2.117773e-03	0.018844548	0.874	0.126
B04	0.0086819802	-1.265781e-03	0.019197647	0.866	0.134
D05	0.0037777625	-6.379054e-03	0.011495393	0.822	0.178
B06	0.0095293505	-4.760244e-03	0.023419874	0.812	0.188
D09	0.0044482653	-4.066690e-03	0.012324899	0.808	0.192
S03	0.0154206703	-1.668538e-02	0.041133011	0.800	0.200
D08	0.0003742537	-4.028928e-04	0.001116957	0.691	0.309
D01	0.0007972038	-4.891676e-03	0.006278674	0.668	0.332

# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()

trait	mean	p1	p9	positive_support	negative_support
B10	-0.0006304998	-0.0009540545	-0.0003081384	0.020	0.980
D02	-0.0151032927	-0.0273709242	-0.0026929023	0.061	0.939
S01	-0.0134498482	-0.0281931644	0.0021095325	0.132	0.868
B07	-0.0093150264	-0.0210303996	0.0025766178	0.142	0.858
S02	-0.0039566793	-0.0121727363	0.0061940272	0.170	0.830
B02	-0.0035889597	-0.0161455535	0.0073469073	0.362	0.638
B08	-0.0017867040	-0.0149536701	0.0101457248	0.533	0.467

positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
#      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")